The life-saving benefit of dexamethasone in severe COVID-19 is linked to a reversal of monocyte dysregulation
Rainer Knoll
07/25/2023
1. R requirements
1.1 Install and load packages
Load the following packages needed for analysis:
packages <- c("reshape2", "stringr", "Seurat", "tibble", "ggplot2", "dplyr", "harmony", "clusterProfiler", "ComplexHeatmap", "ggpubr", "ggrepel", "ggnewscale", "Matrix", "useful", "RColorBrewer", "org.Hs.eg.db", "tidyr", "circlize", "patchwork", "pheatmap", "gridExtra", "grid", "cowplot", "viridis", "fmsb")
lapply(packages, require, character.only = TRUE)1.2 Load functions
Functions for downstream analysis
source("Functions.R")hallmark_genes_symbols <- read.gmt("./files/h.all.v7.4.symbols.gmt")
hallmark_genes <- read.gmt("./files/h.all.v7.4.entrez.gmt")
reactome_genes_symbols <- read.gmt("./files/c2.cp.reactome.v7.4.symbols.gmt")
reactome_genes <- read.gmt("./files/c2.cp.reactome.v7.4.entrez.gmt")
KEGG_genes_symbols<-read.gmt("./files/c2.cp.kegg.v7.4.symbols.gmt")
KEGG_genes<-read.gmt("./files/c2.cp.kegg.v7.4.entrez.gmt")
gmtfile_hallmarks <- hallmark_genes1.3 Load colors
Used colors in this manuscript.
color_clusters<-c(brewer.pal(n = 9,name = "Set1"),
brewer.pal(n = 8,name = "Set2"),
brewer.pal(n = 12,name = "Set3"),
brewer.pal(n = 12,name = "Paired"),
brewer.pal(n = 9,name = "Pastel1"),
brewer.pal(n = 8,name = "Pastel2"),
brewer.pal(n = 8,name = "Accent"))
colors_celltype <- c("CD8 T cells" = "#B02457",
"CD4 T cells" = '#ea7800',
"NK cells" = '#66357D',
"Megakaryocytes" = '#5a6b22',
"B cells" = "#fcc113",
"Monocytes" = "#9dcfcf",
"Plasmablasts" = '#17315f',
"mDCs" = "#C19BA8",
"pDCs" = "#BF40A8",
"Erythrocytes" = "red",
"Prol. cells" = '#9eceaa')
colors_monocytes <- c( "IFIhi" = "#6d3b8f",
"IL1Bhi" = "#b50e3b",
"S100Ahi" = "#4D39D1",
"CXCLhi" = "#fcc113",
"Dexa_response" = "#41b6c4",
"classical" = "#FBB99C",
"THBDhi" = "#e04b07",
"PF4+" = "#e54585",
"HBhi" = "lightgrey",
"non_classical" = "#FF99FF")
colors_severity <- c("severe" = "#CD5C5C",
"moderate" = "#5B9BD5")
colors_severity_ctrl <- c("severe_ctrl" = "indianred",
"moderate_ctrl" = "orange")2. Load Datasets
2.1 PBMC data
seurat_PBMC <- readRDS("COVID_Dexa_seurat_PBMC.RDS")
seurat_PBMC## An object of class Seurat
## 25498 features across 114181 samples within 2 assays
## Active assay: RNA (25486 features, 2000 variable features)
## 1 other assay present: HTO
## 2 dimensional reductions calculated: pca, umap2.2 monocyte data
seurat_monocytes <- readRDS("COVID_Dexa_seurat_monocytes.RDS")
seurat_monocytes## An object of class Seurat
## 25498 features across 23648 samples within 2 assays
## Active assay: RNA (25486 features, 1000 variable features)
## 1 other assay present: HTO
## 3 dimensional reductions calculated: pca, umap, harmony3. Figure 1
3.1 UMAP by celltype
Visualization of cell types on the UMAP embedding.
seurat_PBMC$cluster_labels_res.0.4 <- factor(seurat_PBMC$cluster_labels_res.0.4, levels = c( "Monocytes", "mDCs", "pDCs",
"Megakaryocytes", "CD4 T cells",
"CD8 T cells", "NK cells", "B cells",
"Plasmablasts", "Prol. cells", "Erythrocytes"))
DimPlot(seurat_PBMC, reduction = "umap", group.by = "cluster_labels_res.0.4", cols = colors_celltype, raster = F)
3.2 UMAP by severity
Subsetting the data to controls only (untreated COVID-19) and highlighting severity.
seurat_ctrls <- subset(seurat_PBMC, group %in% c("moderate_ctrl", "severe_ctrl"))
p1 <-DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(14000), sizes.highlight = 0.1, cols.highlight = "#FFB651", cols = "lightgrey")+NoLegend()+ggtitle("moderate")
p2 <- DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(row.names(seurat_ctrls@meta.data[seurat_ctrls@meta.data$group == "severe_ctrl",]), 14000), sizes.highlight = 0.1, cols.highlight = "#CD5C5C", cols = "lightgrey")+NoLegend()+ggtitle("severe")
p1 + p2
3.3 UMAP by early/late
Subsetting the data to controls only (untreated COVID-19) and highlighting timepoints.
seurat_ctrls$timepoints <- ifelse(seurat_ctrls$SO.to.sample > 10, "late", "early" )
p1 <- DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(row.names(seurat_ctrls@meta.data[seurat_ctrls@meta.data$timepoints == "early",]), 17000), sizes.highlight = 0.1, cols.highlight = "#c51b7d", cols = "lightgrey", order = F, raster = F)+NoLegend()+ggtitle("Early <= 10d")
p2 <-DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(row.names(seurat_ctrls@meta.data[seurat_ctrls@meta.data$timepoints == "late",]), 17000), sizes.highlight = 0.1, cols.highlight = "#92c5de", cols = "lightgrey", order = F, raster = F)+NoLegend()+ggtitle("Late >10d")
p1+p2
3.4 Marker by cell type
Top marker by cell types.
seurat_PBMC$cluster_labels_res.0.4 <- factor(seurat_PBMC$cluster_labels_res.0.4, levels = c( "Monocytes", "mDCs", "pDCs",
"Megakaryocytes", "CD4 T cells",
"CD8 T cells", "NK cells", "B cells",
"Plasmablasts", "Prol. cells", "Erythrocytes"))
Idents(seurat_PBMC) <- "cluster_labels_res.0.4"
cluster.markers.wilcox_cluster_labels_res.0.4 <- FindAllMarkers(object = seurat_PBMC,
only.pos = TRUE,
min.pct = 0.2,
logfc.threshold = 0.5,
min.diff.pct = 0.2,
test.use = "wilcox", max.cells.per.ident = 7000, verbose = F
)
top <- cluster.markers.wilcox_cluster_labels_res.0.4 %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DotPlot(seurat_PBMC, features = rev(unique(top$gene)),
group.by = "cluster_labels_res.0.4", col.max = 1.5, col.min = -1.5)+
scale_color_gradientn(colors = rev(RColorBrewer::brewer.pal(n =20, name = "RdBu")))+
theme(axis.text.x = element_text(size = 14, angle = 90, hjust = 1, vjust =0.5),
legend.position = "bottom", axis.text.y = element_text(face = "italic", size = 12))+
coord_flip()
3.5 Monocyte HLA and alarmin
Checking “ground truth” for HLA and alarmin expression in untreated COVID-19 patients monocytes.
seurat_ctrls_monos <- subset(seurat_monocytes, group %in% c("severe_ctrl", "moderate_ctrl"))
# selected HLA
p1 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("HLA-DRA"), pt.size = 0, cols = colors_severity_ctrl)+ stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+theme( axis.title.y = element_blank(), axis.title.x = element_blank())
p2 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("HLA-DRB1"), pt.size = 0, cols = colors_severity_ctrl)+ stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme( axis.title.y = element_blank(), axis.title.x = element_blank())
p3 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("HLA-DPA1"), pt.size = 0, cols = colors_severity_ctrl)+ stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())
# selected alarmins
p4 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("S100A8"), pt.size = 0, cols = colors_severity_ctrl)+ stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p5 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("S100A9"), pt.size = 0, cols = colors_severity_ctrl)+ stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p6 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("S100A12"), pt.size = 0, cols = colors_severity_ctrl)+ stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p1+p2+p3+p4+p5+p6+plot_layout(ncol = 6)
Testing significance for selected genes
df <- FindMarkers(seurat_ctrls_monos, group.by = "group", ident.1 = "severe_ctrl", ident.2 = "moderate_ctrl", features = c("HLA-DRA", "HLA-DRB1", "HLA-DPA1", "S100A8", "S100A9", "S100A12"))
df$gene <- row.names(df)
df[sort(df$gene), c("avg_log2FC", "p_val_adj")]## avg_log2FC p_val_adj
## HLA-DPA1 -0.9061705 6.252979e-208
## HLA-DRA -0.9458845 2.937397e-278
## HLA-DRB1 -1.3566956 0.000000e+00
## S100A12 1.1022834 3.743919e-301
## S100A8 1.0036270 0.000000e+00
## S100A9 0.5324372 2.161001e-1593.6 IFN module monocytes
IFN-gene expression enrichment in monocytes by timepoints.
seurat <- subset(seurat_monocytes, group_severity %in% c("moderate", "severe"))
seurat$order_group <- paste(seurat$order, seurat$group, sep = "_")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_severe_Dexa", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_severe_ctrl", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_moderate_Dexa", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_moderate_ctrl", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "late_severe_Dexa", replacement = "late_Dexa")
seurat$order_group <- gsub(seurat$order_group, pattern = "late_moderate_Dexa", replacement = "late_Dexa")
seurat$order_group <- gsub(seurat$order_group, pattern = "late_severe_ctrl", replacement = "late_Ctrl")
seurat$order_group <- gsub(seurat$order_group, pattern = "late_moderate_ctrl", replacement = "late_Ctrl")
seurat$sev_order_group <- paste(seurat$group_severity, seurat$order_group, sep = "_")
IFNA_signature <- c("IFI6", "ISG15", "IFITM1", "ISG20",
"IFI27", "IFI30", "IFIH1", "IFIT1",
"IFIT2", "IFIT3", "IFITM2", "IFITM3",
"XAF1", "MX1", "MX2")
seurat <- AddModuleScore(seurat, features = list(IFNA_signature), name = "IFNA_signature")
df <- seurat@meta.data[,c("donorID", "sampleID", "sev_order_group", "IFNA_signature1")]
df %>% group_by(sev_order_group) %>% summarize(mean = mean(IFNA_signature1)) -> df.mean
df <- merge(df.mean, df, by = "sev_order_group")
p <- ggplot(df, aes(x = sev_order_group, y = IFNA_signature1))+
geom_violin(aes(fill = mean))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle(label = "IFNA_signature", subtitle = paste("n genes = ", length(IFNA_signature)))+
NoLegend()
p
3.7 PBMC DE call
DEG call for Dexa vs. control for both moderate and severe COVID-19.
seurat <- subset(seurat_PBMC, order == "late" & group_severity %in% c("moderate", "severe"))
# define comparisons
comparisons<-data.frame(comparison=c("sev_Dex_vs_ctrl","mod_Dex_vs_ctrl"),
group1=c("severe_Dexa","moderate_Dexa"),
group2=c("severe_ctrl","moderate_ctrl"))
#define cell types to include
meta<-seurat@meta.data
celltype_to_show <- c("Monocytes", "CD4 T cells", "CD8 T cells", "NK cells", "B cells")
DEgenes.all<-data.frame()
for(i in comparisons$comparison){
groups_oi<-as.vector(comparisons[comparisons$comparison==i,])
seurat_tmp<-subset(seurat, subset = group %in% groups_oi)
for(j in celltype_to_show){
Idents(seurat_tmp)<-"cluster_labels_res.0.4"
seurat_celltype<-subset(seurat_tmp, subset = cluster_labels_res.0.4==j)
Idents(seurat_celltype)<-"group"
tmp <-FindMarkers(object = seurat_celltype,
ident.1 = comparisons[comparisons$comparison==i,]$group1,
ident.2 = comparisons[comparisons$comparison==i,]$group2,
only.pos = FALSE,
min.pct = 0.1,
logfc.threshold = 0.25,
test.use = "wilcox")
#add information
tmp$celltype<-j
tmp$cluster<-paste0(j,"_",i)
tmp$comparison<-i
tmp$celltype_comparison<-paste0(j,"_",i)
tmp$gene<-row.names(tmp)
tmp$regulation<-ifelse(tmp$avg_log2FC>0,"up", "down")
tmp$significance<-ifelse(tmp$p_val_adj<0.05,"<0.05", "ns")
DEgenes.all<-rbind(DEgenes.all,tmp)
}
}
rm(seurat_tmp, seurat_celltype,groups_oi, tmp)DEgenes<-DEgenes.all[DEgenes.all$p_val_adj<0.05,]
DEgenes$celltype <- factor(DEgenes$celltype, levels = c("Monocytes","B cells", "CD4 T cells", "CD8 T cells", "NK cells"))
p <- ggplot(DEgenes, aes(comparison, fill=regulation))+
geom_bar(position = position_dodge(), color = "black")+
scale_fill_manual(values = c(up = "#F59F9F", down = "#76B2ED"))+
scale_alpha_manual(values = c("<0.05"=1, "ns"=0.8))+
theme_linedraw()+
geom_text(stat='count',
aes(label=..count..),
position=position_dodge2(width=0.8), size=3, vjust = -0.3, hjust = 0.6)+
theme(axis.text.x = element_text(angle = 90, hjust=1,vjust = 0))+
ggtitle("DE genes (p.adj<0.05)")+
xlab("")+
ylim(0,260)+
facet_wrap(.~celltype,nrow = 1)+
theme(strip.background =element_rect(fill="grey"))+
theme(strip.text = element_text(color = "black",size = 10))
p
3.8 Common DEG severe
Identification of genes that are DEG in at least 3 cell types in severe COVID-19 Dexa vs. control.
DEgenes_sev <- DEgenes[DEgenes$comparison == "sev_Dex_vs_ctrl", ]
# check commond DEG down
DEgenes_down <- DEgenes_sev[DEgenes_sev$regulation == "down", ]
tmp_all_down <- c()
for(i in unique(DEgenes_down$celltype)){
tmp <- DEgenes_down[DEgenes_down$celltype == i, ]
name <- gsub(pattern = "_sev_Dex_vs_ctrl_down", replacement = "", x = tmp$cluster_regulation)
name <- unique(name)
name <- gsub(name, pattern = " ", replacement = "_")
tmp <- tmp$gene
tmp_all_down <- c(tmp_all_down, tmp)
}
df <- as.data.frame(tmp_all_down)
colnames(df) <- "genes"
genes_multiple_celltypes_down <- vector()
for(i in unique(df$genes)){
tmp <- df[df$genes == i, ]
tmp <- as.data.frame(tmp)
if(nrow(tmp) >= 3){
genes_multiple_celltypes_down <- c(i, genes_multiple_celltypes_down)
}
}
# check common DEG down
DEgenes_sev <- DEgenes[DEgenes$comparison == "sev_Dex_vs_ctrl", ]
DEgenes_up <- DEgenes_sev[DEgenes_sev$regulation == "up", ]
tmp_all_up <- c()
for(i in unique(DEgenes_up$celltype)){
tmp <- DEgenes_up[DEgenes_up$celltype == i, ]
name <- gsub(pattern = "_sev_Dex_vs_ctrl_up", replacement = "", x = tmp$cluster_regulation)
name <- unique(name)
name <- gsub(name, pattern = " ", replacement = "_")
tmp <- tmp$gene
tmp_all_up <- c(tmp_all_up, tmp)
}
df <- as.data.frame(tmp_all_up)
colnames(df) <- "genes"
genes_multiple_celltypes_up <- vector()
for(i in unique(df$genes)){
tmp <- df[df$genes == i, ]
tmp <- as.data.frame(tmp)
if(nrow(tmp) >= 3){
genes_multiple_celltypes_up <- c(i, genes_multiple_celltypes_up)
}
}
list(up = genes_multiple_celltypes_up,
down = genes_multiple_celltypes_down)## $up
## [1] "HLA-F" "SF3A2" "RPS5" "RPL39" "TXNIP" "TSC22D3"
##
## $down
## [1] "ZNF331" "NR4A2" "EEF1A1" "SNORD3A" "EIF5A"
## [6] "AL590867.2" "HMGB2" "RPL9P9" "IFITM1" "MT-ATP8"
## [11] "RPL10" "RPL6" "PTMA" "RPL4" "RPL14"
## [16] "MTND2P28" "LITAF" "AP005202.1" "AC109326.1" "AL365357.1"
## [21] "FTH1"common DEG severe include TSC22D3 and TXNIP (up) and FTH1 and IFITM1 (down).
seurat <- subset(seurat_PBMC, order == "late" & group_severity %in% c("severe"))
seurat <- subset(seurat, cluster_labels_res.0.4 %in% c("Monocytes","B cells", "CD4 T cells",
"CD8 T cells", "NK cells"))
seurat$cluster_labels_res.0.4 <- factor(seurat$cluster_labels_res.0.4 , levels = c("Monocytes","B cells", "CD4 T cells",
"CD8 T cells", "NK cells"))
p1 <- VlnPlot(seurat, features = "TSC22D3", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p2 <- VlnPlot(seurat, features = "TXNIP", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p3 <- VlnPlot(seurat, features = "FTH1", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p4 <- VlnPlot(seurat, features = "IFITM1", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p1+p2+p3+p4+plot_layout(ncol = 4)
statistics check- to identidy DEG.
genes <- c("TSC22D3", "TXNIP", "FTH1", "IFITM1")
tmp <- DEgenes_sev[DEgenes_sev$gene %in% genes, ]
tmp <- tmp[order(tmp$gene, decreasing = F), ]
tmp[, c("celltype", "gene", "p_val_adj")]## celltype gene p_val_adj
## FTH1 Monocytes FTH1 6.317354e-109
## FTH11 CD4 T cells FTH1 5.448513e-34
## FTH12 CD8 T cells FTH1 4.043850e-07
## FTH13 NK cells FTH1 2.339069e-13
## FTH14 B cells FTH1 4.965424e-37
## IFITM1 CD4 T cells IFITM1 8.021860e-35
## IFITM11 CD8 T cells IFITM1 1.947397e-07
## IFITM12 NK cells IFITM1 2.496873e-10
## TSC22D3 Monocytes TSC22D3 5.211175e-130
## TSC22D31 CD4 T cells TSC22D3 5.756550e-55
## TSC22D32 CD8 T cells TSC22D3 2.721599e-06
## TSC22D33 NK cells TSC22D3 6.003117e-03
## TSC22D34 B cells TSC22D3 4.622287e-39
## TXNIP Monocytes TXNIP 2.363646e-41
## TXNIP1 CD8 T cells TXNIP 1.668797e-09
## TXNIP2 NK cells TXNIP 3.626118e-05
## TXNIP3 B cells TXNIP 2.128502e-423.9 GOEA PBMC
3.9.1 Common GOEA in multiple cell types
DEgenes$cluster_regulation<-paste(DEgenes$celltype_comparison, DEgenes$regulation, sep="_")
DE.list<-list()
for(i in unique(DEgenes$cluster_regulation)){
DE.list[[i]]<-DEgenes[DEgenes$cluster_regulation==i,]$gene
}
GOEA_celltypes_oi<-GOEA.of.list(seurat_object= seurat,
DE.list,
GeneSets =c("GO","Hallmark", "KEGG", "Reactome"
),
GOntology = "BP",
pCorrection = "bonferroni", # choose the p-value adjustment method (Steffi: bonferroni)
pvalueCutoff = 0.05, # set the unadj. or adj. p-value cutoff (depending on correction method)
qvalueCutoff = 0.05 # set the q-value cutoff (FDR corrected)
)
GOEA_celltypes_oi$KEGGresults <- NULL
GOEA_celltypes_oi$HALLMARKresults$database <- "Hallmark"
GOEA_celltypes_oi$GOresults$database <- "GO"
GOEA_celltypes_oi$REACTOMEresults$database <- "Reactome"
GOEA_celltypes_oi$REACTOMEresults$Description <- gsub(GOEA_celltypes_oi$REACTOMEresults$Description, pattern = "REACTOME_", replacement = "")
GOEA_all <- rbind(GOEA_celltypes_oi$GOresults,
rbind(GOEA_celltypes_oi$REACTOMEresults, GOEA_celltypes_oi$HALLMARKresults))
GOEA_all$geneSet <- row.names(GOEA_all)
GOEA_all$group <- sapply(strsplit(GOEA_all$geneSet, split = "_"), `[`, 2)
GOEA_all %>% group_by(Condition) %>% arrange(p.adjust, .by_group = TRUE) -> GOEA_all
GOEA_all %>% group_by(Condition) %>% arrange(Count, .by_group = TRUE) %>% top_n(n= 10, wt = p.adjust) %>% pull(Description) ->terms
GOEA_all <- GOEA_all[GOEA_all$Description %in% terms,]
GOEA_all$Description <- ifelse(nchar(GOEA_all$Description) > 50,
paste(substr(GOEA_all$Description, 1, 50), "[...]", sep=""),
GOEA_all$Description)
GOEA_all$Description <- factor(GOEA_all$Description, levels = rev(unique(GOEA_all$Description)))
#change color code of dots
GOEA_all$DE.regulation<-"no"
if(nrow(GOEA_all[grep("_up", GOEA_all$Condition),])>0){
GOEA_all[grep("_up", GOEA_all$Condition),]$DE.regulation<-"up"
}
if(nrow(GOEA_all[grep("_down", GOEA_all$Condition),])>0){
GOEA_all[grep("_down", GOEA_all$Condition),]$DE.regulation<-"down"
}
GOEA_all$Cond_short <- gsub(GOEA_all$Condition, pattern = "_.*", replacement = "")
GOEA_all_small <- GOEA_all[, c("Cond_short", "Description", "Count", "group", "DE.regulation")]
nrow(GOEA_all_small)## [1] 302df <- GOEA_all_small %>% group_by(Description, group, DE.regulation) %>% summarize(count = n())
GOEA_all_small <- GOEA_all_small[GOEA_all_small$Description %in% df[df$count >= 3, ]$Description, ]
GOEA_all_small <- GOEA_all_small[GOEA_all_small$Cond_short %in% c("Monocytes", "B cells", "CD4 T cells", "CD8 T cells", "NK cells"), ]
GOEA_all_small$Cond_short <- factor(GOEA_all_small$Cond_short, levels = c("Monocytes", "B cells", "CD4 T cells", "CD8 T cells", "NK cells"))
p <- ggplot(data = GOEA_all_small ,aes(x=Cond_short, y=Description, size=Count)) +
geom_point(fill = "black",pch=21)+theme_classic()+
theme_linedraw()+
theme(axis.text.y = element_text(color = "black"),
axis.text.x = element_text(angle=90, vjust=0.5,hjust=1, color = "black"),
text = element_text(size = 12))+
scale_fill_manual(values= c(sev = "#CD5C5C", mod = "#FFB651"))+
xlab("")+ylab("")+
scale_radius(breaks = waiver())+
facet_wrap(DE.regulation ~ group, scales = "fixed")
p
3.9.2 NFKB down spider
selection of HALLMARK_TNFA_SIGNALING_VIA_NFKB as the best term for down regulated genes and checking enrichment again in a spyder plot.
df <- GOEA_all_small[GOEA_all_small$Description %in% "HALLMARK_TNFA_SIGNALING_VIA_NFKB" & GOEA_all_small$DE.regulation == "down", ]
extension <- data.frame("Cond_short" = c("NK cells", "CD4 T cells", "CD8 T cells",
"CD4 T cells", "NK cells"),
"Description" = c("HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"HALLMARK_TNFA_SIGNALING_VIA_NFKB", "HALLMARK_TNFA_SIGNALING_VIA_NFKB"),
"Count" = c(0, 0, 0, 0, 0),
"group" = c("mod", "mod", "mod", "sev", "sev"),
"DE.regulation" = c("down", "down", "down", "down", "down"))
df <- rbind(df, extension)
df$Cond_short <- factor(df$Cond_short, levels = c("Monocytes", "B cells", "CD4 T cells", "CD8 T cells", "NK cells"))
df_mod <- df[df$group == "mod", ]
df_sev <- df[df$group == "sev", ]
max <- ceiling(max(c(max(df_sev$Count), max(df_mod$Count))))
data <- matrix(rep(0, 5), ncol = 5)
colnames(data) <- c("Monocytes", "B cells", "NK cells", "CD8 T cells", "CD4 T cells")
tmp <- rbind(rep(40, 5), data, df_mod$Count[match(colnames(data), df_mod$Cond_short)])
tmp <- as.data.frame(tmp)
tmp[is.na(tmp)] <- 0
radarchart(tmp, axistype=1 ,
#custom polygon
pcol=rgb(1, 0.714, 0.318,0.9) , pfcol=rgb(1, 0.714, 0.318,0.5) , plwd=4 ,
#custom the grid
cglcol="darkgrey", cglty=1, axislabcol="black", cglwd=1, caxislabels=seq(0,40,10), seg = 4, title = "moderate NF-kB TF prediciton",
#custom labels
vlcex=1 )
data <- matrix(rep(0, 5), ncol = 5)
colnames(data) <- c("Monocytes", "B cells", "NK cells", "CD8 T cells", "CD4 T cells")
tmp <- rbind(rep(40, 5), data, df_sev$Count[match(colnames(data), df_sev$Cond_short)])
tmp <- as.data.frame(tmp)
tmp[is.na(tmp)] <- 0
radarchart(tmp, axistype=1 ,
#custom polygon
pcol=rgb(0.804, 0.361, 0.361, 0.9) , pfcol=rgb(0.804, 0.361, 0.361,0.5) , plwd=4 ,
#custom the grid
cglcol="darkgrey", cglty=1, axislabcol="black", cglwd=1, caxislabels=seq(0,40,10), seg = 4, title = "severate NF-kB TF prediciton",
#custom labels
vlcex=1 )
4. Figure 2
4.1 FC-FC plot Dexa core
To identify the Core Dexa signature in monocytes we compare DEG from the comparison Dexa vs control for moderate and severe COVID-19.
DEG_list <- list()
for(i in c("moderate", "severe")){
tmp <- subset(seurat_monocytes, subset = group_severity == i)
for (x in c("late")) {
tmp_sample_order <- subset(tmp, subset = order == x)
tmp_sample_order <- SetIdent(tmp_sample_order, value = "group")
ctrl <- unique(tmp_sample_order$group)[grep(unique(tmp_sample_order$group), pattern = "ctrl")]
idx <- ifelse(grep(unique(tmp_sample_order$group), pattern = "ctrl") == 1, 2, 1)
treat <- unique(tmp_sample_order$group)[idx]
tryCatch({
DEG<-FindMarkers(tmp_sample_order,
ident.1 = treat,
ident.2 = ctrl,
min.pct = 0.1,
logfc.threshold = 0,
test.use = "wilcox",
slot = "data",
only.pos = F,
verbose = F)
DEG$Gene <- rownames(DEG)
DEG$group_severity <- paste(i)
DEG$sample_comparison <- paste(x)
DEG$comparison <- paste0(i, "_", x, "_treatment_vs_ctrl")
DEG_list[[paste0(i, "_", x, "_treatment_vs_ctrl")]]<-DEG
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}
}
tmp1 <- DEG_list$moderate_late_treatment_vs_ctrl[intersect(row.names(DEG_list$moderate_late_treatment_vs_ctrl), row.names(DEG_list$severe_late_treatment_vs_ctrl)),]
tmp2 <- DEG_list$severe_late_treatment_vs_ctrl[intersect(row.names(DEG_list$moderate_late_treatment_vs_ctrl), row.names(DEG_list$severe_late_treatment_vs_ctrl)),]
colnames(tmp1) <- paste(colnames(tmp1), "moderate", sep = "_")
colnames(tmp2) <- paste(colnames(tmp2), "severe", sep = "_")
data <- cbind(tmp1, tmp2)
data$Gene <- row.names(data)
data_rank_dexa_mod_severe <- data
data_rank_dexa_mod_severe$regulation <- ifelse(data_rank_dexa_mod_severe$p_val_adj_moderate < 0.05 & data_rank_dexa_mod_severe$p_val_adj_severe < 0.05, "DE_both", "DE_one")
data_rank_dexa_mod_severe$keep <- ifelse(data_rank_dexa_mod_severe$p_val_adj_moderate >= 0.05 & data_rank_dexa_mod_severe$p_val_adj_severe >= 0.05, "remove", "keep")
data_rank_dexa_mod_severe[data_rank_dexa_mod_severe$keep %in% "remove", "regulation"] <- "n.s."
df <- data_rank_dexa_mod_severe
down_label <- df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_both",]
down_label_mod <- down_label[order(down_label$avg_log2FC_moderate), "Gene"][1:15]
down_label_sev <- down_label[order(down_label$avg_log2FC_severe), "Gene"][1:15]
down_label_gene <- unique(down_label_sev, down_label_mod)
up_label <- df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_both",]
up_label_mod <- up_label[order(up_label$avg_log2FC_moderate, decreasing = T), "Gene"][1:15]
up_label_sev <- up_label[order(up_label$avg_log2FC_severe, decreasing = T), "Gene"][1:15]
up_label_gene <- unique(up_label_sev, up_label_mod)
p <- ggplot(df, aes(x = avg_log2FC_moderate, y = avg_log2FC_severe))+
geom_point(data =df,
pch = 21,size = 2, color = "black", fill = "white")+
geom_point(data= df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_one",],
aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#c6dbef")+
geom_point(data= df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_one",],
aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#fcc5c0")+
geom_text_repel(data= df[df$Gene %in% down_label_gene, ],
aes(x=avg_log2FC_moderate, y=avg_log2FC_severe, label = Gene),
size = 4, color = "black", fontface = 'italic', max.overlaps = 40)+
geom_point(data= df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_both",],
aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#2171b5")+
geom_text_repel(data= df[df$Gene %in% up_label_gene, ],
aes(x=avg_log2FC_moderate, y=avg_log2FC_severe, label = Gene),
size = 4, color = "black", fontface = 'italic', max.overlaps = 40)+
geom_point(data= df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_both",], aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#cb181d")+
theme_bw()+
theme_linedraw()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 10),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 10),
strip.background = element_rect(fill=c("lightgrey")),
strip.text = element_text(colour = "black", face = "bold", size = 12),
plot.title = element_text(face = "bold", size = 12))+
xlab("log2FC moderate")+
ylab("log2FC severe")+
ggtitle("monocytes dexa: severe vs moderate")+
geom_hline(yintercept = 0)+
geom_vline(xintercept = 0)+
geom_abline(linetype = 2)+
xlim(c(-1.6,3.7))+
ylim(c(-1.7, 2.1))
p
4.2 Dexa-induced gene expression violins
seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat_monocytes_late, group_severity %in% c("moderate", "severe"))
seurat$group <- factor(seurat$group, levels = c("moderate_ctrl", "moderate_Dexa",
"severe_ctrl", "severe_Dexa"))
p1<-VlnPlot(seurat, features = c("TSC22D3"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)
p2<-VlnPlot(seurat, features = c("IL1R2"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p3<-VlnPlot(seurat, features = c("CD163"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p4<-VlnPlot(seurat, features = c("SAP30"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p5<- VlnPlot(seurat, features = c("PER1"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p6<- VlnPlot(seurat, features = c("ZFP36L2"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p1+p2+p3+p4+p5+p6+plot_layout(ncol = 3)
checking significance for DEG:
tmp <- DEG_list$moderate_late_treatment_vs_ctrl
tmp <- tmp[tmp$Gene %in% c("TSC22D3", "IL1R2", "CD163", "SAP30", "PER1", "ZFP36L2"), c("group_severity", "Gene", "p_val_adj", "avg_log2FC")]
tmp[order(tmp$Gene), ]## group_severity Gene p_val_adj avg_log2FC
## CD163 moderate CD163 5.191492e-122 1.1818631
## IL1R2 moderate IL1R2 0.000000e+00 2.6483134
## PER1 moderate PER1 1.735679e-78 0.8029548
## SAP30 moderate SAP30 4.544743e-135 1.2053391
## TSC22D3 moderate TSC22D3 3.557426e-91 0.7122654
## ZFP36L2 moderate ZFP36L2 7.352414e-68 0.9031032tmp <- DEG_list$severe_late_treatment_vs_ctrl
tmp <- tmp[tmp$Gene %in% c("TSC22D3", "IL1R2", "CD163", "SAP30", "PER1", "ZFP36L2"), c("group_severity", "Gene", "p_val_adj", "avg_log2FC")]
tmp[order(tmp$Gene), ]## group_severity Gene p_val_adj avg_log2FC
## CD163 severe CD163 2.929259e-20 0.5129026
## IL1R2 severe IL1R2 5.249462e-87 1.7486982
## PER1 severe PER1 1.472531e-07 0.3160355
## SAP30 severe SAP30 5.188564e-23 0.4672367
## TSC22D3 severe TSC22D3 2.212959e-117 0.6846452
## ZFP36L2 severe ZFP36L2 2.804138e-08 0.45059094.3 Core Heatmap
Core Dexa DEG as heatmap (top 20 up and down).
df <- data_rank_dexa_mod_severe
down_label <- df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_both",]
down_label_mod <- down_label[order(down_label$avg_log2FC_moderate), "Gene"][1:20]
down_label_sev <- down_label[order(down_label$avg_log2FC_severe), "Gene"][1:20]
down_label_gene <- unique(down_label_sev, down_label_mod)
up_label <- df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_both",]
up_label_mod <- up_label[order(up_label$avg_log2FC_moderate, decreasing = T), "Gene"][1:20]
up_label_sev <- up_label[order(up_label$avg_log2FC_severe, decreasing = T), "Gene"][1:20]
up_label_gene <- unique(up_label_sev, up_label_mod)
df_1 <- DEG_list$moderate_late_treatment_vs_ctrl[row.names(DEG_list$moderate_late_treatment_vs_ctrl) %in% c(up_label_gene, down_label_gene),]
colnames(df_1) <- paste(colnames(df_1), "moderate", sep = "_")
df_2 <- DEG_list$severe_late_treatment_vs_ctrl[row.names(DEG_list$severe_late_treatment_vs_ctrl) %in% c(up_label_gene, down_label_gene),]
colnames(df_2) <- paste(colnames(df_2), "severe", sep = "_")
df_1 <- df_1[c(up_label_gene, down_label_gene),]
df_2 <- df_2[c(up_label_gene, down_label_gene),]
df_all <- cbind(df_1, df_2)
df_all$Gene <- row.names(df_all)
df_melt <- melt(df_all[,c("Gene", "avg_log2FC_moderate", "avg_log2FC_severe")])
df_melt$variable <- gsub(pattern = "avg_log2FC_", replacement = "", df_melt$variable)
df_melt$Gene <- factor(df_melt$Gene, levels = c(down_label_gene, up_label_gene))
df_melt$variable <- factor(df_melt$variable, levels = c("moderate", "severe"))
tmp <- df_all[, c("avg_log2FC_moderate", "avg_log2FC_severe")]
ord <- hclust( dist(tmp, method = "euclidean"), method = "ward.D" )$order
tmp$gene <- row.names(tmp)
tmp <- melt(tmp[rev(ord), ])
tmp$variable <- gsub(pattern = "avg_log2FC_", replacement = "", tmp$variable)
tmp$gene <- factor(tmp$gene, levels = rev(unique(tmp$gene)))
heatmap <- ggplot(tmp, aes(x = variable, y = gene, fill = value))+
geom_tile(color = "black")+
theme_classic() +
scale_fill_gradient2(midpoint = 0, low = "#2171b5", high = "#cb181d")+
theme(text = element_text(size = 14, color = "black"),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(size = 16, angle = 60, hjust = 1),
plot.title = element_text(size = 16),
strip.text.y = element_text(angle = 0))+
ggtitle("Avg logFC")
heatmap
4.4 Functional enrichment of core genes
plot_list <- list()
down <- df[df$avg_log2FC_moderate <(-0.1) & df$avg_log2FC_severe <(-0.1) & df$regulation %in% "DE_both", "Gene"]
up <- df[df$avg_log2FC_moderate >0.1 & df$avg_log2FC_severe >0.1 & df$regulation %in% "DE_both", "Gene"]
GSEA_core <- GSEA2(object = seurat_monocytes_late,
up_reg_genes = up, down_reg_genes = down,
name = "Dexa vs. Ctrl.",
GeneSets =c("GO", "Hallmark", "KEGG"),
GOntology = "BP",
pval.use = "p_val_adj",
pCorrection = "bonferroni", # choose the p-value adjustment method
pvalueCutoff = 0.05,
qvalueCutoff = 0.05 # set the q-value cutoff (FDR corrected)
)
GSEA_core$GOup$database <- "GO"
GSEA_core$HALLMARKup$database <- "HALLMARK"
df_up <- rbind(GSEA_core$GOup,
GSEA_core$HALLMARKup)
df_up <- df_up[order(df_up$Count, decreasing = T),]
df_up$Description <- factor(df_up$Description, levels = rev(unique(df_up$Description)))
df_up$Description <- gsub(df_up$Description, pattern = "HALLMARK_", replacement = "")
df_up <- df_up[order(df_up$Count, decreasing = T),]
df_up$Description <- factor(df_up$Description, levels = rev(unique(df_up$Description)))
GSEA_core$GOdown$database <- "GO"
GSEA_core$HALLMARKdown$database <- "HALLMARK"
GSEA_core$GOdown$comparison <- NULL
df_down <- rbind(GSEA_core$GOdown,
GSEA_core$HALLMARKdown)
df_down <- df_down[order(df_down$Count, decreasing = T),]
df_down$Description <- factor(df_down$Description, levels = rev(unique(df_down$Description)))
df_down$Description <- gsub(df_down$Description, pattern = "HALLMARK_", replacement = "")
df_down <- df_down[order(df_down$Count, decreasing = T),]
df_down$Description <- factor(df_down$Description, levels = rev(unique(df_down$Description)))
df_up$comparison <- NULL
df_up$reg <- "up"
df_down$reg <- "down"
df_both <- rbind(df_up[df_up$Count >0, ], head(df_down[order(df_down$Count, decreasing = T), ], 5))
df_both$reg <- factor(df_both$reg, levels = c("up", "down"))
df_both$Description <- factor(df_both$Description, levels = rev(unique(df_both$Description)) )
p <- ggplot(df_both, aes(x = reg, y = Description, fill = reg)) +
geom_point(aes(size = Count), pch = 21, color = "black")+
ggtitle(label = "functional Dexa core", subtitle = "(p<0.05)")+
theme_classic()+
scale_fill_manual(values = c("#cb181d", "#2171b5"))+
theme_bw()+
theme_linedraw()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12),
legend.text = element_text( size = 12),
plot.title = element_text(face = "bold", size = 14))+
xlab("")+
ylab("Description")
p
4.5 Wang GC signature monocytes
enrichment of GC in-vitro signatures from Wang et al.
seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat_monocytes_late, group_severity %in% c("moderate", "severe"))
GC_wang <- read.csv("./files/Wang_et_al_Nature_GC_treatment_monocytes.csv", sep = ";")
GC_wang <- na.omit(GC_wang)
GC_wang_up <- GC_wang[order(GC_wang$Mo.DMSO_vs_TA_log2FoldChange, decreasing = T),][1:100,]
GC_wang_up <- GC_wang_up$gene
GC_wang_down <- GC_wang[order(GC_wang$Mo.DMSO_vs_TA_log2FoldChange, decreasing = F),][1:100,]
GC_wang_down <- GC_wang_down$gene
seurat <- AddModuleScore(object = seurat,
features = list(GC_wang_up),
name = 'GC_wang_up')
seurat <- AddModuleScore(object = seurat,
features = list(GC_wang_down),
name = 'GC_wang_down')
df <- seurat@meta.data[,c("donorID", "sampleID", "group", "GC_wang_up1")]
df %>% group_by(group) %>% summarize(mean = mean(GC_wang_up1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "group")
df <- df[order(df$scale_mean, decreasing = T),]
df$group <- factor(df$group, levels = c("moderate_ctrl", "moderate_Dexa", "severe_ctrl", "severe_Dexa"))
p1 <- ggplot(df, aes(x = group, y = GC_wang_up1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("GC_wang_up")+
NoLegend()
df <- seurat@meta.data[,c("donorID", "sampleID", "group", "GC_wang_down1")]
df %>% group_by(group) %>% summarize(mean = mean(GC_wang_down1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "group")
df <- df[order(df$scale_mean, decreasing = T),]
df$group <- factor(df$group, levels = c("moderate_ctrl", "moderate_Dexa", "severe_ctrl", "severe_Dexa"))
p2 <- ggplot(df, aes(x = group, y = GC_wang_down1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("GC_wang_down")+
NoLegend()
p1+p2
4.6 Monocyte UMAP
monocyte_levels <- c("IFIhi", "IL1Bhi", "S100Ahi", "CXCLhi", "Dexa_response", "classical", "THBDhi", "PF4+", "HBhi", "non_classical")
seurat_monocytes$cluster_labels_res.0.7 <- factor(seurat_monocytes$cluster_labels_res.0.7, levels = monocyte_levels)
DimPlot(seurat_monocytes, reduction = "umap", cols = colors_monocytes, group.by = "cluster_labels_res.0.7", order = F)
4.7 Monocyte state marker plot
seurat_monocytes$cluster_labels_res.0.7 <- factor(seurat_monocytes$cluster_labels_res.0.7, levels = monocyte_levels)
Idents(seurat_monocytes) <- "cluster_labels_res.0.7"
cluster.markers.wilcox_mono <- FindAllMarkers(object = seurat_monocytes,
only.pos = TRUE,
min.pct = 0.2,
logfc.threshold = 0.25,
test.use = "wilcox", verbose = F
)
top <- cluster.markers.wilcox_mono %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
genes <- c("IFI6", "IFI27", "IFI44L", "XAF1", "MX1",
"IL1B", "CCL3", "CCL3L3", "CCL4", "CCL4L2", "TNF",
"S100A12", "S100A8", "S100A9", "RETN", "LGALS1",
"CXCL1", "CXCL2", "CXCL3", "CXCL8", "MAFB", "HBEGF",
"AREG", "IL1R2", "TSC22D3", "THBS1", "SAP30", "JDP2", "ZFP36L2", "EREG", "CXCR4", "CD163", "FKBP5", "KLF9", "PER1",
"HLA-DPA1", "HLA-DOB1", "HLA-DRA", "HLA-DRB1", "HLA-DQB1",
"RGCC", "LMNA", "PHLDA1", "THBD", "CCL7", "EMP1",
"PPBP", "PF4", "TUBB1", "NRGN", "SPARC",
"HBB", "HBA1", "HBA2",
"FCGR3A", "CX3CR1", "CDKN1C", "C1QA", "C1QB", "C1QC")
seurat_monocytes$cluster_labels_res.0.7 <- factor(seurat_monocytes$cluster_labels_res.0.7, levels = rev(monocyte_levels))
p <- DotPlot(seurat_monocytes,
group.by = "cluster_labels_res.0.7",
features = genes)+
RotatedAxis()+
scale_color_gradientn(colors = rev(RColorBrewer::brewer.pal(n =20, name = "RdBu")))+
theme(axis.title = element_blank()) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p
4.8 Comparison to Schulte-Schrepping marker list
p <- DotPlot(seurat_monocytes,
group.by = "cluster_labels_res.0.7",
features = c("MTND2P28", "HLA-DQB1", "HLA-DRB1", "HLA-DMB", "HLA-DPB1", "HLA-DRB5", # Classical Monocytes
"S100A12", "S100A8", "S100A9", "CXCL8", "ALOX5AP", "IRS2", "PLBD1", # HLA-DRlo S100Ahi
"NRGN", "PF4", "PBPP", "MYL9", "CAVIN2", "SPARC", "TUBB1", "ITGA2B", # HLA-DRlo PF4+
"IFI6", "IFI27", "MX1", "ISG15", "SELL", "CD163", "OAS2", "RETN", # HLA-DRlo CD163hi
"HBB", "HBA2", "HBA1", "HBD", # HLA-DRhi HBBhi
"RGCC", "PHLDA1", "EMP1", "THBD", "SKG1", "ATP2B1", "CD83", # HLA-DRhi CD83hi
"FCGR3A", "CDKN1C", "C1QB", "CX3CR1", "MS4A7", "C1QA", "C1QC" ))+ # non-classical Monocytes
RotatedAxis()+
scale_color_gradientn(colors = rev(RColorBrewer::brewer.pal(n =20, name = "RdBu")))+
theme(axis.title = element_blank()) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, face = "italic"))
p
4.9 Signature enrichment by monocyte state
response_glucocorticoid <- read.delim("./files/GO_term_summary_20230303_050304.txt", row.names = NULL)
response_glucocorticoid <- response_glucocorticoid$MGI.Gene.Marker.ID
response_glucocorticoid <- toupper(response_glucocorticoid)
response_steroid_hormone <- clusterProfiler::read.gmt("./files/GO_response_to_steroid_hormone.gmt")
response_steroid_hormone <- response_steroid_hormone$gene
seurat <- seurat_monocytes
seurat <- AddModuleScore(object = seurat,
features = list(response_glucocorticoid),
name = 'response_glucocorticoid')
seurat <- AddModuleScore(object = seurat,
features = list(response_steroid_hormone),
name = 'response_steroid')
seurat <- AddModuleScore(object = seurat,
features = list(GC_wang_up),
name = 'GC_wang_up')
seurat <- AddModuleScore(object = seurat,
features = list(GC_wang_down),
name = 'GC_wang_down')
df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "GC_wang_up1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(GC_wang_up1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")
df <- df[order(df$scale_mean, decreasing = T),]
levels_plot <- unique(df$cluster_labels_res.0.7)
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)
p1 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = GC_wang_up1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("GC_wang_up")+
NoLegend()
df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "GC_wang_down1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(GC_wang_down1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")
df <- df[order(df$scale_mean, decreasing = T),]
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)
p2 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = GC_wang_down1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("GC_wang_down")+
NoLegend()
df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "response_glucocorticoid1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(response_glucocorticoid1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)
p3 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = response_glucocorticoid1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("response_glucocorticoid")+
NoLegend()
df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "response_steroid1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(response_steroid1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)
p4 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = response_steroid1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 14, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("response_steroid")+
NoLegend()
p1+p2+p3+p4+plot_layout(ncol = 2)
5. Figure 3
5.1 Dexa response state distribution
seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat_quant <- subset(seurat_monocytes_late, group_severity %in% c("moderate", "severe"))
# quantify cells of each sample per cluster
cells_cluster <- as.matrix(confusionMatrix(paste0(seurat_quant$cluster_labels_res.0.7),
paste0(seurat_quant$sampleID)))
percentages <- round(t(t(cells_cluster)/colSums(cells_cluster))*100,3)
# transpose data for inspection
meta <- data.frame(seurat_quant@meta.data)
COVID_Dexa.percentages <- melt(percentages)
colnames(COVID_Dexa.percentages) <- c("cluster_labels_res.0.7","sampleID","percent")
idx <- match(COVID_Dexa.percentages$sampleID,meta$sampleID)
COVID_Dexa.percentages$group <- meta$group[idx]
COVID_Dexa.percentages$outcome <- meta$outcome[idx]
data <- COVID_Dexa.percentages[COVID_Dexa.percentages$cluster_labels_res.0.7 == "Dexa_response",]
data$group <- factor(data$group, levels = c("moderate_ctrl", "moderate_Dexa",
"severe_ctrl", "severe_Dexa"))
p <- ggplot(data, aes(fill = group, y=percent, x=group)) +
geom_boxplot(outlier.size = 0, outlier.color = "white")+
geom_point(pch = 21, color = "black", fill = "gray30" ,size=2,position = position_jitter(width = 0.1, height = 0))+
theme_classic()+
ggtitle("Dexa response state")+
scale_fill_manual(values= c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61")) +
xlab("")+
theme(axis.text.x = element_text(angle = 45, hjust =1, size = 12),
axis.text.y = element_text(size = 12),
plot.title = element_text(size=16, face = "bold"))+
ggpubr::stat_compare_means(comparisons = list(c("severe_ctrl", "severe_Dexa"),
c("moderate_ctrl", "moderate_Dexa")), method = "wilcox.test")
p
5.2 Dexa response state by outcome
seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat_quant <- subset(seurat_monocytes_late, group %in% c("severe_Dexa"))
# quantify cells of each sample per cluster
cells_cluster <- as.matrix(confusionMatrix(paste0(seurat_quant$cluster_labels_res.0.7),
paste0(seurat_quant$sampleID)))
percentages <- round(t(t(cells_cluster)/colSums(cells_cluster))*100,3)
# transpose data for inspection
meta <- data.frame(seurat_quant@meta.data)
COVID_Dexa.percentages <- melt(percentages)
colnames(COVID_Dexa.percentages) <- c("cluster_labels_res.0.7","sampleID","percent")
idx <- match(COVID_Dexa.percentages$sampleID,meta$sampleID)
COVID_Dexa.percentages$group <- meta$group[idx]
COVID_Dexa.percentages$outcome <- meta$outcome[idx]
data <- COVID_Dexa.percentages[COVID_Dexa.percentages$cluster_labels_res.0.7 == "Dexa_response",]
data$outcome <- factor(data$outcome, levels = c("deceased", "survived"))
p <- ggplot(data, aes(fill = outcome, y=percent, x=outcome)) +
geom_boxplot(outlier.size = 0, outlier.color = "white")+
geom_point(pch = 21, color = "black", fill = "gray30" ,size=2,position = position_jitter(width = 0.1, height = 0))+
theme_classic()+
ggtitle("Dexa response state")+
scale_fill_manual(values= c("#5B7674", "#57B391")) +
xlab("")+
theme(axis.text.x = element_text(angle = 45, hjust =1, size = 12),
axis.text.y = element_text(size = 12),
plot.title = element_text(size=16, face = "bold"))+
ggpubr::stat_compare_means(comparisons = list(c("deceased", "survived")), method = "wilcox.test")
p
## 5.3 Volcano severe survivors Dexa vs control
DEG_list <- list()
for(i in c("severe")){
tmp <- subset(seurat_monocytes, subset = group_severity == i)
tmp <- subset(tmp, subset = outcome == "survived")
for (x in c("late")) {
tmp_sample_order <- subset(tmp, subset = order == x)
tmp_sample_order <- SetIdent(tmp_sample_order, value = "group")
ctrl <- unique(tmp_sample_order$group)[grep(unique(tmp_sample_order$group), pattern = "ctrl")]
idx <- ifelse(grep(unique(tmp_sample_order$group), pattern = "ctrl") == 1, 2, 1)
treat <- unique(tmp_sample_order$group)[idx]
tryCatch({
DEG<-FindMarkers(tmp_sample_order,
ident.1 = treat,
ident.2 = ctrl,
min.pct = 0.1,
logfc.threshold = 0,
test.use = "wilcox",
slot = "data",
only.pos = F,
verbose = F)
DEG$Gene <- rownames(DEG)
DEG$group_severity <- paste(i)
DEG$sample_comparison <- paste(x)
DEG$comparison <- paste0(i, "_", x, "_treatment_vs_ctrl")
DEG_list[[paste0(i, "_", x, "_treatment_vs_ctrl")]]<-DEG
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}
}
up_reg <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC>0,]
up_reg <- up_reg[order(up_reg$avg_log2FC, -up_reg$p_val_adj, decreasing = T),]$Gene
down_reg <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$avg_log2FC, -down_reg$p_val_adj, decreasing = F),]$Gene
up_reg_pval <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC>0,]
up_reg_pval <- up_reg_pval[order(up_reg_pval$p_val_adj, decreasing = F),]$Gene
down_reg <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$p_val_adj, decreasing = F),]$Gene
p <-ggplot(DEG_list$severe_late_treatment_vs_ctrl,aes(x=avg_log2FC,y=-log10(p_val_adj))) +
geom_point(pch = 21, fill = "darkgrey", color = "black", size = 2) +
geom_text_repel(data=DEG_list$severe_late_treatment_vs_ctrl[c("CST3", "CD74", "HLA-DRB1", "HLA-DRA",
"AREG", "HLA-DPA1", "TSC22D3",
"THBS1", "CEBPD", "HLA-DPB1",
"IL1R2", "MT-ND3", "HLA-F", "CXCR4",
"ZFP36L2"),],
max.overlaps = 20,
aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),
size=3, color = "black")+
geom_text_repel(data=DEG_list$severe_late_treatment_vs_ctrl[c("CCL4","S100A9", "S100A8", "S100A12",
"CCL3", "CCL4L2", "CCL3L3", "IL1B",
"NAMPTP1","IL1RN", "RETN", "CLU",
"CXCL3", "G0S2", "CXCL2", "PHLDA1",
"RGCC", "CXCL8", "HP", "BCL2A1", "CD36",
"FTH1", "CST7"),],
max.overlaps = 20,
aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),size=3,
color = "black")+
geom_point(data=DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC>0.2,],
aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#57B391", size = 2, pch = 21, color = "black")+
geom_point(data=DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC<c(-0.2),],
aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#F08FB6", size = 2, pch = 21, color = "black")+
ggtitle(paste0("survivor DEG analysis (severe COVID-19)"))+
theme_classic()+
theme_bw()+
theme_linedraw()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 8),
axis.text.y = element_text( size = 8),
axis.title.x = element_text( size = 8),
axis.text.x = element_text( size = 8),
legend.title = element_blank(),
legend.position = "none",
plot.title = element_text(face = "bold", size = 10))+
xlab("log2FC")+
ylab("-log10(p_val_adj)")+
geom_hline(yintercept = 0)+
geom_vline(xintercept = 0)+
xlim(c(-2, 2.5))+
ylim(c(0, 220))
p 
5.4 volcano of selected DEG Dexa survivors
seurat <- subset(seurat_monocytes, order == "late" & outcome == "survived" & group_severity %in% c("severe"))
p1 <- VlnPlot(seurat, features = "HLA-DRB1", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p2 <- VlnPlot(seurat, features = "HLA-DRA", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p3 <- VlnPlot(seurat, features = "HLA-DPA1", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p4 <- VlnPlot(seurat, features = "CD74", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p5 <- VlnPlot(seurat, features = "S100A8", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p6 <- VlnPlot(seurat, features = "S100A9", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p7 <- VlnPlot(seurat, features = "S100A12", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p8 <- VlnPlot(seurat, features = "CCL4", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p9 <- VlnPlot(seurat, features = "IL1B", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p10 <- VlnPlot(seurat, features = "CXCL3", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p11 <- VlnPlot(seurat, features = "CXCL2", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p12 <- VlnPlot(seurat, features = "CXCL8", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
NoLegend()+
theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p1+p2 + p5+p6 + p9+p10 +
p3+p4 + p7+p8 + p11+p12 + plot_layout(nrow = 2, ncol = 6)
5.5 S100A state distribution
seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat_quant <- subset(seurat_monocytes_late, group %in% c("severe_Dexa"))
# quantify cells of each sample per cluster
cells_cluster <- as.matrix(confusionMatrix(paste0(seurat_quant$cluster_labels_res.0.7),
paste0(seurat_quant$sampleID)))
percentages <- round(t(t(cells_cluster)/colSums(cells_cluster))*100,3)
# transpose data for inspection
meta <- data.frame(seurat_quant@meta.data)
COVID_Dexa.percentages <- melt(percentages)
colnames(COVID_Dexa.percentages) <- c("cluster_labels_res.0.7","sampleID","percent")
idx <- match(COVID_Dexa.percentages$sampleID,meta$sampleID)
COVID_Dexa.percentages$group <- meta$group[idx]
COVID_Dexa.percentages$outcome <- meta$outcome[idx]
data <- COVID_Dexa.percentages[COVID_Dexa.percentages$cluster_labels_res.0.7 == "S100Ahi",]
data$outcome <- factor(data$outcome, levels = c("deceased", "survived"))
p <- ggplot(data, aes(fill = outcome, y=percent, x=outcome)) +
geom_boxplot(outlier.size = 0, outlier.color = "white")+
geom_point(pch = 21, color = "black", fill = "gray30" ,size=2,position = position_jitter(width = 0.1, height = 0))+
theme_classic()+
ggtitle("S100Ahi state")+
scale_fill_manual(values= c("#5B7674", "#57B391")) +
xlab("")+
theme(axis.text.x = element_text(angle = 45, hjust =1, size = 12),
axis.text.y = element_text(size = 12),
plot.title = element_text(size=16, face = "bold"))+
ggpubr::stat_compare_means(comparisons = list(c("deceased", "survived")), method = "wilcox.test")
p
5.6 Dexa outcome DEG volcano (severe COVID-19)
seurat <- subset(seurat_monocytes, subset = group_severity == "severe" & order == "late")
seurat <- SetIdent(seurat, value = "group_outcome")
DEG_deceased <- list()
x= "late"
tryCatch({
DEG<-FindMarkers(seurat,
ident.1 = "severe_Dexa_survived",
ident.2 = "severe_Dexa_deceased",
min.pct = 0.1,
logfc.threshold = 0,
test.use = "wilcox",
slot = "data",
only.pos = F,
verbose = F)
DEG$Gene <- rownames(DEG)
DEG_deceased[[ paste0(x, "_treatment_alive_vs_deceased")]]<-DEG
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
df <- DEG_deceased$late_treatment_alive_vs_deceased
up_reg <- df[df$p_val_adj<.05&df$avg_log2FC>0,]
up_reg <- up_reg[order(up_reg$avg_log2FC, -up_reg$p_val_adj, decreasing = T),]$Gene
down_reg <- df[df$p_val_adj<.05&df$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$avg_log2FC, -down_reg$p_val_adj, decreasing = F),]$Gene
up_reg_pval <- df[df$p_val_adj<.05&df$avg_log2FC>0,]
up_reg_pval <- up_reg_pval[order(up_reg_pval$p_val_adj, decreasing = F),]$Gene
down_reg <- df[df$p_val_adj<.05&df$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$p_val_adj, decreasing = F),]$Gene
p <- ggplot(df,aes(x=avg_log2FC,y=-log10(p_val_adj))) +
geom_point(pch = 21, fill = "darkgrey", color = "black", size = 2) +
geom_text_repel(data=df[c("S100A8", "S100A9", "S100A12", "PKM", "CLU",
"LGALS1", "MCEMP1", "RETN", "MT-ATP8",
"TMEM176B", "PLAC8", "SERPINB2"),],
max.overlaps = 20,
aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),size=3, color = "black")+
geom_text_repel(data=df[c("HLA-DRB1", "CD74", "HLA-DRA",
"HLA-DPA1", "CEBPD", "HLA-DPB1",
"ZFP36L2", "HLA-DQB1", "IFI6",
"AREG", "CCL3","CCL4", "CCL3L3",
"IL1R2", "HLA-DMB", "HLA-DRB5",
"RPL34", "CST3", "RPS13"),],
max.overlaps = 20,
aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),size=3, color = "black")+
geom_point(data=df[df$p_val_adj<.05&df$avg_log2FC>0.2,],
aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#57B391", size = 2, pch = 21, color = "black")+
geom_point(data=df[df$p_val_adj<.05&df$avg_log2FC<c(-0.2),],
aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#5B7674", size = 2, pch = 21, color = "black")+
ggtitle(paste0("Dexa outcome DEG analysis (severe COVID-19)"))+
theme_classic()+
theme_bw()+
theme_linedraw()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 8),
axis.text.y = element_text( size = 8),
axis.title.x = element_text( size = 8),
axis.text.x = element_text( size = 8),
legend.title = element_blank(),
legend.position = "none",
plot.title = element_text(face = "bold", size = 10))+
xlab("log2FC")+
ylab("-log10(p_val_adj)")+
geom_hline(yintercept = 0)+
geom_vline(xintercept = 0)+
xlim(c(-1.3, 1.6))+
ylim(c(0, 280))
p
5.7 selected outcome DEG
seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat_monocytes_late, group_severity %in% "severe")
seurat <- subset(seurat, group_outcome %in% "severe_ctrl_deceased", invert = T)
seurat$group_outcome <- factor(seurat$group_outcome, levels = rev(c("severe_ctrl_survived",
"severe_Dexa_deceased",
"severe_Dexa_survived")))
p <- DotPlot(seurat, group.by = "group_outcome",
features = c("S100A8", "S100A9", "S100A12",
"CD74", "HLA-DRB1", "HLA-DRA",
"HLA-DPA1", "HLA-DPB1", "HLA-DQB1",
"HLA-DMB", "ZFP36L2", "IFI6", "AREG"),
col.max = 1.5, col.min = -1.5)+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
scale_color_gradient2(low = "blue", mid = "white", high = "red")
p
check for correlation of the selected genes.
seurat$group_outcome_donor <- paste(seurat$group_outcome, seurat$donorID, sep = "_")
df <- AverageExpression(seurat, features = c("S100A8", "S100A9", "S100A12",
"CD74", "HLA-DRB1", "HLA-DRA", "HLA-DPA1", "HLA-DPB1", "HLA-DQB1", "HLA-DMB", "ZFP36L2", "IFI6", "AREG"
), group.by = "group_outcome")$RNA
df## severe_Dexa_survived severe_Dexa_deceased severe_ctrl_survived
## S100A8 119.475170 206.7166969 203.6020996
## S100A9 103.012413 158.3877037 172.6169736
## S100A12 11.206291 22.2545702 17.8782846
## CD74 45.761604 21.9769300 27.6624037
## HLA-DRB1 12.116669 3.6801748 5.7532979
## HLA-DRA 17.651813 7.1160404 10.2649664
## HLA-DPA1 4.711761 1.3969531 2.4362415
## HLA-DPB1 2.601729 0.7854829 1.2535002
## HLA-DQB1 1.841001 0.5625530 1.5580807
## HLA-DMB 1.578108 0.6232262 0.8807122
## ZFP36L2 4.901420 2.1062185 2.2731095
## IFI6 6.430200 2.1234939 3.8872294
## AREG 4.219793 1.2213803 2.0480760df <- as.data.frame(t(scale(t(df))))
df <- as.data.frame(df)
cor_df <- cor(df)
cor_df <- cor_df[c("severe_ctrl_survived", "severe_Dexa_deceased", "severe_Dexa_survived"), ]
pheatmap(cor_df[,"severe_ctrl_survived"],
cluster_cols = F,
cluster_rows = F,
display_numbers = T,
number_color = "black")
5.8 HLA-DRlo S100Ahi Schulte-Schrepping signature enrichment
seurat_Schulte <- readRDS("./files/seurat_COVID19_monocytes_jonas_FG_2020-07-29.rds")
seurat_Schulte@meta.data$cell_ID <- row.names(seurat_Schulte@meta.data)
seurat_Schulte$dataset <- "Schulte"
seurat_Schulte$cluster_labels_res.0.2<- seurat_Schulte$RNA_snn_res.0.2
Idents(seurat_Schulte) <- "cluster_labels_res.0.2"
seurat_Schulte <- RenameIdents(object = seurat_Schulte,
`0` = "HLA-DRlo_S100Ahi",
`1` = "HLA-DRhi_CD83hi",
`2` = "HLA-DRlo_CD163hi",
`3` = "Classical_Monocytes",
`4` = "Non-classical_Monocytes",
`5` = "HLA-DRhi_HBBhi",
`6` = "HLA-DRlo_PF4+")
seurat_Schulte@meta.data$cluster_labels_res.0.2 <- Idents(seurat_Schulte)
marker <- FindMarkers(seurat_Schulte, ident.1 = "HLA-DRlo_S100Ahi",
group.by = "cluster_labels_res.0.2",
logfc.threshold = 0.25, min.pct = 0.2, only.pos = T)
marker <- marker[marker$p_val_adj < 0.05, ]
marker <- marker[order(marker$avg_log2FC, decreasing = T), ]
Schulte_S100_top50 <- head(row.names(marker), 50)
seurat <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat, group_severity %in% c("severe"))
seurat <- AddModuleScore(object = seurat,
features = list(Schulte_S100_top50),
name = 'HLA_DRlo_S100Ahi_signature')
df <- seurat@meta.data[,c("donorID", "sampleID", "group_outcome", "HLA_DRlo_S100Ahi_signature1")]
df %>% group_by(group_outcome) %>% summarize(mean = mean(HLA_DRlo_S100Ahi_signature1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "group_outcome")
df <- df[order(df$scale_mean, decreasing = T),]
df$group_outcome <- factor(df$group_outcome, levels = c("severe_ctrl_deceased",
"severe_ctrl_survived",
"severe_Dexa_deceased",
"severe_Dexa_survived"))
p <- ggplot(df, aes(x = group_outcome, y = HLA_DRlo_S100Ahi_signature1))+
geom_violin(aes(fill = scale_mean))+
stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = -0.5)+
theme_classic()+
theme(panel.grid.major = element_line("white"),
panel.grid.minor = element_line("white"),
axis.title.y = element_text( size = 12),
axis.text.y = element_text( size = 12),
axis.title.x = element_text( size = 12),
axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
legend.text = element_text( size = 12),
plot.title = element_text(size = 10, hjust = 0.5))+
xlab("")+
ylab("Module Score")+
ggtitle("HLA-DRlo S100Ahi signature (Schulte-Schrepping, n=50)")
p
6. Enviornment info
R.version## _
## platform x86_64-pc-linux-gnu
## arch x86_64
## os linux-gnu
## system x86_64, linux-gnu
## status
## major 4
## minor 1.0
## year 2021
## month 05
## day 18
## svn rev 80317
## language R
## version.string R version 4.1.0 (2021-05-18)
## nickname Camp PontanezensessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 grid stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] fmsb_0.7.5 viridis_0.6.1 viridisLite_0.4.0
## [4] cowplot_1.1.1 gridExtra_2.3 pheatmap_1.0.12
## [7] patchwork_1.1.1 circlize_0.4.13 tidyr_1.1.4
## [10] org.Hs.eg.db_3.13.0 AnnotationDbi_1.54.1 IRanges_2.26.0
## [13] S4Vectors_0.30.2 Biobase_2.52.0 BiocGenerics_0.38.0
## [16] RColorBrewer_1.1-2 useful_1.2.6 Matrix_1.5-4.1
## [19] ggnewscale_0.4.8 ggrepel_0.9.1 ggpubr_0.4.0
## [22] ComplexHeatmap_2.8.0 clusterProfiler_4.0.5 harmony_0.1.0
## [25] Rcpp_1.0.7 dplyr_1.0.7 ggplot2_3.3.5
## [28] tibble_3.1.5 SeuratObject_4.1.3 Seurat_4.3.0
## [31] stringr_1.4.0 reshape2_1.4.4
##
## loaded via a namespace (and not attached):
## [1] scattermore_0.7 bit64_4.0.5 knitr_1.33
## [4] irlba_2.3.3 data.table_1.14.2 KEGGREST_1.32.0
## [7] RCurl_1.98-1.5 doParallel_1.0.16 generics_0.1.0
## [10] RSQLite_2.2.8 shadowtext_0.0.9 RANN_2.6.1
## [13] future_1.12.0 bit_4.0.4 enrichplot_1.12.3
## [16] spatstat.data_3.0-1 httpuv_1.6.3 assertthat_0.2.1
## [19] xfun_0.24 hms_1.1.1 jquerylib_0.1.4
## [22] evaluate_0.14 promises_1.2.0.1 fansi_0.5.0
## [25] readxl_1.3.1 igraph_1.2.6 DBI_1.1.1
## [28] htmlwidgets_1.5.4 spatstat.geom_3.2-1 purrr_0.3.4
## [31] ellipsis_0.3.2 backports_1.2.1 deldir_1.0-5
## [34] vctrs_0.3.8 Cairo_1.5-12.2 ROCR_1.0-11
## [37] abind_1.4-5 cachem_1.0.6 withr_2.4.2
## [40] ggforce_0.3.3 progressr_0.13.0 sctransform_0.3.5
## [43] treeio_1.16.2 goftest_1.2-3 cluster_2.1.2
## [46] DOSE_3.18.3 ape_5.5 lazyeval_0.2.2
## [49] crayon_1.4.1 spatstat.explore_3.2-1 pkgconfig_2.0.3
## [52] labeling_0.4.2 tweenr_1.0.2 GenomeInfoDb_1.28.4
## [55] vipor_0.4.5 nlme_3.1-152 rlang_1.1.1
## [58] globals_0.14.0 lifecycle_1.0.1 miniUI_0.1.1.1
## [61] downloader_0.4 ggrastr_0.2.3 cellranger_1.1.0
## [64] polyclip_1.10-0 matrixStats_0.61.0 lmtest_0.9-38
## [67] aplot_0.1.1 carData_3.0-4 zoo_1.8-9
## [70] beeswarm_0.4.0 ggridges_0.5.3 GlobalOptions_0.1.2
## [73] png_0.1-7 rjson_0.2.20 bitops_1.0-7
## [76] KernSmooth_2.23-20 Biostrings_2.60.2 blob_1.2.2
## [79] shape_1.4.6 qvalue_2.24.0 spatstat.random_3.1-5
## [82] rstatix_0.7.0 gridGraphics_0.5-1 ggsignif_0.6.3
## [85] scales_1.1.1 memoise_2.0.0 magrittr_2.0.1
## [88] plyr_1.8.6 ica_1.0-2 zlibbioc_1.38.0
## [91] compiler_4.1.0 scatterpie_0.1.7 clue_0.3-60
## [94] fitdistrplus_1.1-6 cli_3.6.1 XVector_0.32.0
## [97] listenv_0.8.0 pbapply_1.5-0 MASS_7.3-54
## [100] tidyselect_1.1.1 stringi_1.7.5 forcats_0.5.1
## [103] highr_0.9 yaml_2.2.1 GOSemSim_2.18.1
## [106] sass_0.4.0 fastmatch_1.1-3 tools_4.1.0
## [109] future.apply_1.2.0 rio_0.5.27 rstudioapi_0.13
## [112] foreach_1.5.1 foreign_0.8-81 farver_2.1.0
## [115] Rtsne_0.15 ggraph_2.0.5 digest_0.6.28
## [118] shiny_1.7.1 car_3.0-11 broom_0.7.9
## [121] later_1.3.0 RcppAnnoy_0.0.19 httr_1.4.2
## [124] colorspace_2.0-2 tensor_1.5 reticulate_1.22
## [127] splines_4.1.0 uwot_0.1.14 yulab.utils_0.0.4
## [130] tidytree_0.3.5 spatstat.utils_3.0-3 graphlayouts_0.7.1
## [133] sp_1.6-0 ggplotify_0.1.0 plotly_4.10.0
## [136] xtable_1.8-4 jsonlite_1.7.2 ggtree_3.0.4
## [139] tidygraph_1.2.0 ggfun_0.0.4 R6_2.5.1
## [142] pillar_1.6.3 htmltools_0.5.2 mime_0.12
## [145] glue_1.4.2 fastmap_1.1.0 BiocParallel_1.26.2
## [148] codetools_0.2-18 fgsea_1.18.0 utf8_1.2.2
## [151] lattice_0.20-44 bslib_0.3.1 spatstat.sparse_3.0-1
## [154] ggbeeswarm_0.6.0 curl_4.3.2 leiden_0.3.9
## [157] zip_2.2.0 GO.db_3.13.0 openxlsx_4.2.4
## [160] survival_3.2-11 limma_3.48.3 rmarkdown_2.11
## [163] munsell_0.5.0 DO.db_2.9 GetoptLong_1.0.5
## [166] GenomeInfoDbData_1.2.6 iterators_1.0.13 haven_2.4.3
## [169] gtable_0.3.0